home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / C / Applications / Python 1.3.3 / Python 133 SRC / Modules / cmathmodule.c < prev    next >
Text File  |  1996-01-12  |  5KB  |  315 lines
  1. /* Complex math module */
  2.  
  3. /* much code borrowed from mathmodule.c */
  4.  
  5. #include "allobjects.h"
  6. #include "complexobject.h"
  7.  
  8. #include <errno.h>
  9.  
  10. #include "mymath.h"
  11.  
  12. #ifdef i860
  13. /* Cray APP has bogus definition of HUGE_VAL in <math.h> */
  14. #undef HUGE_VAL
  15. #endif
  16.  
  17. #ifdef HUGE_VAL
  18. #define CHECK(x) if (errno != 0) ; \
  19.     else if (-HUGE_VAL <= (x) && (x) <= HUGE_VAL) ; \
  20.     else errno = ERANGE
  21. #else
  22. #define CHECK(x) /* Don't know how to check */
  23. #endif
  24.  
  25. #ifndef M_PI
  26. #define M_PI (3.141592653589793239)
  27. #endif
  28.  
  29. /* First, the C functions that do the real work */
  30.  
  31. /* constants */
  32. static complex c_1 = {1., 0.};
  33. static complex c_half = {0.5, 0.};
  34. static complex c_i = {0., 1.};
  35. static complex c_i2 = {0., 0.5};
  36. static complex c_mi = {0., -1.};
  37. static complex c_pi2 = {M_PI/2., 0.};
  38.  
  39. /* forward declarations */
  40. complex c_log();
  41. complex c_prodi();
  42. complex c_sqrt();
  43.  
  44.  
  45. complex c_acos(x)
  46.     complex x;
  47. {
  48.     return c_neg(c_prodi(c_log(c_sum(x,c_prod(c_i,
  49.             c_sqrt(c_diff(c_1,c_prod(x,x))))))));
  50. }
  51.  
  52. complex c_acosh(x)
  53.     complex x;
  54. {
  55.     return c_log(c_sum(x,c_prod(c_i,
  56.             c_sqrt(c_diff(c_1,c_prod(x,x))))));
  57. }
  58.  
  59. complex c_asin(x)
  60.     complex x;
  61. {
  62.     return c_neg(c_prodi(c_log(c_sum(c_prod(c_i,x),
  63.             c_sqrt(c_diff(c_1,c_prod(x,x)))))));
  64. }
  65.  
  66. complex c_asinh(x)
  67.     complex x;
  68. {
  69.     return c_neg(c_log(c_diff(c_sqrt(c_sum(c_1,c_prod(x,x))),x)));
  70. }
  71.  
  72. complex c_atan(x)
  73.     complex x;
  74. {
  75.     return c_prod(c_i2,c_log(c_quot(c_sum(c_i,x),c_diff(c_i,x))));
  76. }
  77.  
  78. complex c_atanh(x)
  79.     complex x;
  80. {
  81.     return c_prod(c_half,c_log(c_quot(c_sum(c_1,x),c_diff(c_1,x))));
  82. }
  83.  
  84. complex c_cos(x)
  85.     complex x;
  86. {
  87.     complex r;
  88.     r.real = cos(x.real)*cosh(x.imag);
  89.     r.imag = -sin(x.real)*sinh(x.imag);
  90.     return r;
  91. }
  92.  
  93. complex c_cosh(x)
  94.     complex x;
  95. {
  96.     complex r;
  97.     r.real = cos(x.imag)*cosh(x.real);
  98.     r.imag = sin(x.imag)*sinh(x.real);
  99.     return r;
  100. }
  101.  
  102. complex c_exp(x)
  103.     complex x;
  104. {
  105.     complex r;
  106.     double l = exp(x.real);
  107.     r.real = l*cos(x.imag);
  108.     r.imag = l*sin(x.imag);
  109.     return r;
  110. }
  111.  
  112. complex c_log(x)
  113.     complex x;
  114. {
  115.     complex r;
  116.     double l = hypot(x.real,x.imag);
  117.     r.imag = atan2(x.imag, x.real);
  118.     r.real = log(l);
  119.     return r;
  120. }
  121.  
  122. complex c_log10(x)
  123.     complex x;
  124. {
  125.     complex r;
  126.     double l = hypot(x.real,x.imag);
  127.     r.imag = atan2(x.imag, x.real)/log(10.);
  128.     r.real = log10(l);
  129.     return r;
  130. }
  131.  
  132. complex c_prodi(x)
  133.      complex x;
  134. {
  135.     complex r;
  136.     r.real = -x.imag;
  137.     r.imag = x.real;
  138.     return r;
  139. }
  140.  
  141. complex c_sin(x)
  142.     complex x;
  143. {
  144.     complex r;
  145.     r.real = sin(x.real)*cosh(x.imag);
  146.     r.imag = cos(x.real)*sinh(x.imag);
  147.     return r;
  148. }
  149.  
  150. complex c_sinh(x)
  151.     complex x;
  152. {
  153.     complex r;
  154.     r.real = cos(x.imag)*sinh(x.real);
  155.     r.imag = sin(x.imag)*cosh(x.real);
  156.     return r;
  157. }
  158.  
  159. complex c_sqrt(x)
  160.     complex x;
  161. {
  162.     complex r;
  163.     double s,d;
  164.     if (x.real == 0. && x.imag == 0.)
  165.         r = x;
  166.     else {
  167.         s = sqrt(0.5*(fabs(x.real) + hypot(x.real,x.imag)));
  168.         d = 0.5*x.imag/s;
  169.         if (x.real > 0.) {
  170.             r.real = s;
  171.             r.imag = d;
  172.         }
  173.         else if (x.imag >= 0.) {
  174.             r.real = d;
  175.             r.imag = s;
  176.         }
  177.         else {
  178.             r.real = -d;
  179.             r.imag = -s;
  180.         }
  181.     }
  182.     return r;
  183. }
  184.  
  185. complex c_tan(x)
  186.     complex x;
  187. {
  188.     complex r;
  189.     double sr,cr,shi,chi;
  190.     double rs,is,rc,ic;
  191.     double d;
  192.     sr = sin(x.real);
  193.     cr = cos(x.real);
  194.     shi = sinh(x.imag);
  195.     chi = cosh(x.imag);
  196.     rs = sr*chi;
  197.     is = cr*shi;
  198.     rc = cr*chi;
  199.     ic = -sr*shi;
  200.     d = rc*rc + ic*ic;
  201.     r.real = (rs*rc+is*ic)/d;
  202.     r.imag = (is*rc-rs*ic)/d;
  203.     return r;
  204. }
  205.  
  206. complex c_tanh(x)
  207.     complex x;
  208. {
  209.     complex r;
  210.     double si,ci,shr,chr;
  211.     double rs,is,rc,ic;
  212.     double d;
  213.     si = sin(x.imag);
  214.     ci = cos(x.imag);
  215.     shr = sinh(x.real);
  216.     chr = cosh(x.real);
  217.     rs = ci*shr;
  218.     is = si*chr;
  219.     rc = ci*chr;
  220.     ic = si*shr;
  221.     d = rc*rc + ic*ic;
  222.     r.real = (rs*rc+is*ic)/d;
  223.     r.imag = (is*rc-rs*ic)/d;
  224.     return r;
  225. }
  226.  
  227.  
  228. /* And now the glue to make them available from Python: */
  229.  
  230. static object *
  231. math_error()
  232. {
  233.     if (errno == EDOM)
  234.         err_setstr(ValueError, "math domain error");
  235.     else if (errno == ERANGE)
  236.         err_setstr(OverflowError, "math range error");
  237.     else
  238.         err_errno(ValueError); /* Unexpected math error */
  239.     return NULL;
  240. }
  241.  
  242. static object *
  243. math_1(args, func)
  244.     object *args;
  245.     complex (*func) FPROTO((complex));
  246. {
  247.     complex x;
  248.     if (!PyArg_ParseTuple(args, "D", &x))
  249.         return NULL;
  250.     errno = 0;
  251.     x = (*func)(x);
  252.     CHECK(x.real);
  253.     CHECK(x.imag);
  254.     if (errno != 0)
  255.         return math_error();
  256.     else
  257.         return newcomplexobject(x);
  258. }
  259.  
  260. #define FUNC1(stubname, func) \
  261.     static object * stubname(self, args) object *self, *args; { \
  262.         return math_1(args, func); \
  263.     }
  264.  
  265. FUNC1(cmath_acos, c_acos)
  266. FUNC1(cmath_acosh, c_acosh)
  267. FUNC1(cmath_asin, c_asin)
  268. FUNC1(cmath_asinh, c_asinh)
  269. FUNC1(cmath_atan, c_atan)
  270. FUNC1(cmath_atanh, c_atanh)
  271. FUNC1(cmath_cos, c_cos)
  272. FUNC1(cmath_cosh, c_cosh)
  273. FUNC1(cmath_exp, c_exp)
  274. FUNC1(cmath_log, c_log)
  275. FUNC1(cmath_log10, c_log10)
  276. FUNC1(cmath_sin, c_sin)
  277. FUNC1(cmath_sinh, c_sinh)
  278. FUNC1(cmath_sqrt, c_sqrt)
  279. FUNC1(cmath_tan, c_tan)
  280. FUNC1(cmath_tanh, c_tanh)
  281.  
  282.  
  283. static struct methodlist cmath_methods[] = {
  284.     {"acos", cmath_acos, 1},
  285.     {"acosh", cmath_acosh, 1},
  286.     {"asin", cmath_asin, 1},
  287.     {"asinh", cmath_asinh, 1},
  288.     {"atan", cmath_atan, 1},
  289.     {"atanh", cmath_atanh, 1},
  290.     {"cos", cmath_cos, 1},
  291.     {"cosh", cmath_cosh, 1},
  292.     {"exp", cmath_exp, 1},
  293.     {"log", cmath_log, 1},
  294.     {"log10", cmath_log10, 1},
  295.     {"sin", cmath_sin, 1},
  296.     {"sinh", cmath_sinh, 1},
  297.     {"sqrt", cmath_sqrt, 1},
  298.     {"tan", cmath_tan, 1},
  299.     {"tanh", cmath_tanh, 1},
  300.     {NULL,        NULL}        /* sentinel */
  301. };
  302.  
  303. void
  304. initcmath()
  305. {
  306.     object *m, *d, *v;
  307.     
  308.     m = Py_InitModule("cmath", cmath_methods);
  309.     d = getmoduledict(m);
  310.     dictinsert(d, "pi", v = newfloatobject(atan(1.0) * 4.0));
  311.     DECREF(v);
  312.     dictinsert(d, "e", v = newfloatobject(exp(1.0)));
  313.     DECREF(v);
  314. }
  315.